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Abstract. The purpose of this survey chapter is to present a transformation technique that can 
be used in analysis and numerical computation of the early exercise boundary for an American style 
of vanilla options that can be modelled by class of generalized Black-Scholes equations. We an- 
alyze qualitatively and quantitatively the early exercise boundary for a linear as well as a class of 
nonlinear Black-Scholes equations with a volatility coefficient which can be a nonlinear function 
of the second derivative of the option price itself. A motivation for studying the nonlinear Black- 
Scholes equation with a nonlinear volatility arises from option pricing models taking into account 
e.g. nontrivial transaction costs, investor's preferences, feedback and illiquid markets effects and 
risk from a volatile (unprotected) portfolio. We present a method how to transform the free bound- 
ary problem for the early exercise boundary position into a solution of a time depending nonlinear 
nonlocal parabolic equation defined on a fixed domain. We furthermore propose an iterative numer- 
ical scheme that can be used in order to find an approximation of the free boundary. In the case of 
a linear Black-Scholes equation we are able to derive a nonlinear integral equation for the position 
of the free boundary. We present results of numerical approximation of the early exercise boundary 
for various types of linear and nonlinear Black-Scholes equations and we discuss dependence of the 
free boundary on model parameters. Finally, we discuss an application of the transformation method 
for the pricing equation for American type of Asian options. 
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1 Introduction 

According to the classical theory due to Black, Scholes and Merton the price of an option 
in an idealized financial market can be computed from a solution to the well-known Black- 
Scholes linear parabolic equation derived by Black and Scholes in [7], and, independently 
by Merton (see also Kwok [38], Dewynne et al. [16], Hull [30], Wilmott et al. [53]). Recall 
that a European call (put) option is the right but not obligation to purchase (sell) an under- 
lying asset at the expiration price E at the expiration time T. Assuming that the underlying 
asset S follows a geometric Brownian motion 

dS=(g- q)Sdt + aSdW, (1) 

where g is a drift, q is the asset dividend yield rate, a is the volatility of the asset and W 
is the standard Wiener process (cf. [38]), one can derive a governing partial differential 
equation for the price of an option. We remind ourselves that the equation for option's price 
V(S, t) is the following parabolic PDE: 

where a is the volatility of the underlying asset price process, r > is the interest rate of a 
zero-coupon bond, q > is the dividend yield rate. A solution V = V(S, t) represents the 
price of an option if the price of an underlying asset is S > at time t 6 [0, T]. 

The case when the diffusion coefficient a > in Q is constant represents a classical 
Black-Scholes equation originally derived by Black and Scholes in [7]. On the other hand, 
if we assume the volatility coefficient a > to be a function of the solution V itself then 
with such a diffusion coefficient represents a nonlinear generalization of the Black-Scholes 
equation. It is a purpose of this chapter to focus our attention to the case when the diffusion 
coefficient a 2 may depend on the time T — t to expiry, the asset price S and the second 
derivative d$V of the option price (hereafter referred to as T), i.e. 

a = a{S 2 d 2 s V, S,T-t). (3) 
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A motivation for studying the nonlinear Black-Scholes equation (O with a volatility a 
having a general form (O arises from option pricing models taking into account nontrivial 
transaction costs, market feedbacks and/or risk from a volatile (unprotected) portfolio. Re- 
call that the linear Black-Scholes equation with constant a has been derived under several 
restrictive assumptions like e.g. frictionless, liquid and complete markets, etc. We also 
recall that the linear Black-Scholes equation provides a perfectly replicated hedging port- 
folio. In the last decades some of these assumptions have been relaxed in order to model, 
for instance, the presence of transaction costs (see e.g. Leland [39], Hoggard et al. [29], 
Avellaneda and Paras [4]), feedback and illiquid market effects due to large traders choos- 
ing given stock-trading strategies (Frey and Patie [22], Frey and Stremme [23], During 
et al. [17], Schonbucher and Wilmott [49]), imperfect replication and investor's prefer- 
ences (Barles and Soner [8]), risk from unprotected portfolio (Kratka [37], Jandacka and 
Sevcovic [32] or [47]). One of the first nonlinear models is the so-called Leland model 
(cf. [39]) for pricing call and put options under the presence of transaction costs. It has 
been generalized for more complex option strategies by Hoggard, Whaley and Wilmott 
in [29]. In this model the volatility a is given by 

a 2 (S 2 d 2 s V,S,r) = a 2 (l+Lesgn(d 2 s V)) (4) 

where a > is a constant historical volatility of the underlying asset price process and 
Le > is the so-called Leland constant given by Le = /(ayAt) where C > is 

a constant round trip transaction cost per unit dollar of transaction in the assets market and 
At > is the time-lag between portfolio adjustments. 

Notice that dependence of volatility adjustment on the second derivative of the price 
is quite natural. Indeed, in the idealized Black-Scholes theory, the optimal hedge is equal 
to ±dsV and therefore one may expect more frequent transaction in regions with the high 
second derivative dgV (cf. [7]). 

A popular nonlinear generalization of the Black-Scholes equation has been proposed 
by Avellaneda, Levy, and Paras [5] for description of incomplete markets and uncertain but 
bounded volatility. In their model we have 

•><.#4V,S,t)-{% (5) 

where o\ and a<i represent a lower and upper a-priori bound on the otherwise unspecified 
asset price volatility. 

If transaction costs are taken into account perfect replication of the contingent claim 
is no longer possible and further restrictions are needed in the model. By assuming that 
investor's preferences are characterized by an exponential utility function Barles and Soner 
(cf. [8]) derived a nonlinear Black-Scholes equation with the volatility a given by 

a 2 (S 2 d 2 s V,S,T) = & 2 (1 + ^(a 2 e rT S 2 d 2 s V)) (6) 



where ^ is a solution to the ODE: <f'(x) = (V(x) + l)/(2- v /x^ r (x) - x), *(0) = 0, and 
a > is a given constant representing risk aversion. Notice that &(x) = 0(x3 ) for x — ► 
and ty(x) = 0{x) for x — > oo. 
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Another popular model has been derived for the case when the asset dynamics takes into 
account the presence of feedback and illiquid market effects. Frey and Stremme (cf. [22, 
23]) introduced directly the asset price dynamics in the case when a large trader chooses a 
given stock-trading strategy (see also [49]). The diffusion coefficient a is again nonconstant 
and it can be expressed as: 

a 2 (S 2 d 2 s V,S,r) = a 2 {l- Q\(S)Sd 2 s V)~ 2 (7) 

where a 2 , g > are constants and X(S) is a strictly convex function, X(S) > 1. Interest- 
ingly enough, explicit solutions to the Black-Scholes equation with varying volatility as in 
f7]) have been derived by Bordag and Chankova [9] and Bordag and Frey [10]. 

The last example of the Black-Scholes equation with a nonlinearly depending volatility 
is the so-called Risk Adjusted Pricing Methodology model proposed by Kratka in [37] and 
revisited by Jandacka and Sevcovic in [32]. In order to maintain (imperfect) replication of 
a portfolio by the delta hedge one has to make frequent portfolio adjustments leading to a 
substantial increase in transaction costs. On the other hand, rare portfolio adjustments may 
lead to an increase of the risk arising from a volatile (unprotected) portfolio. In the RAPM 
model the aim is to optimize the time-lag At between consecutive portfolio adjustments. 
By choosing At > in such way that the sum of the rate of transaction costs and the rate of 
a risk from unprotected portfolio is minimal one can find the optimal time lag At > 0. In 
the RAPM model, it turns out that the volatility is again nonconstant and it has the following 
form: 

a 2 (S 2 d 2 s V,S,r) =a 2 (l + fx(Sd 2 s V)^ . (8) 

Here & 2 > is a constant historical volatility of the asset price returns and \x = 
3(C 2 .R/27r)3 where C, R > are nonnegative constant representing the transaction cost 
measure and the risk premium measure, resp. (see [32] for details). 

Notice that all the above mentioned nonlinear models are consistent with the original 
Black-Scholes equation in the case the additional model parameters (e.g. Le, a, g, /i) are 
vanishing. If plain call or put vanilla options are concerned then the function V(S, t) is 
convex in S variable and therefore each of the above mentioned models has a diffusion 
coefficient strictly larger than a 2 leading to a larger values of computed option prices. They 
can be therefore identified with higher Ask option prices, i.e. offers to sell an option. 
Furthermore, these models have been considered and analyzed mostly for European style 
of options, i.e. options that can be exercised only at the maturity t = T. On the other hand, 
American options are much more common in financial markets as they allow for exercising 
of an option anytime before the expiry T. In the case of an American call option a solution 
to equation (f2]) is defined on a time dependent domain < S < St(t), < t < T. It is 
subject to the boundary conditions 

V(0,t) = 0, V(S f (t),t)=S f (t)-E, d s V(S f (t),t) = l, (9) 

and terminal pay-off condition at expiry t = T 



V{S,T) = max(S-£,0) 



(10) 
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where E > is a strike price (cf. [16, 38]). One of important problems in this field is 
the analysis of the early exercise boundary Sf(t) and the optimal stopping time (an inverse 
function to Sf(t)) for American call (or put) options on stocks paying a continuous dividend 
yield with a rate q > (or q > 0). However, an exact analytical expression for the free 
boundary profile is not even known for the case when the volatility a is constant. Many 
authors have investigated various approximation models leading to approximate expressions 
for valuing American call and put options: analytic approximations (Barone-Adesi and 
Whaley [6], Kuske and Keller [36], Dewynne et al. [16], Geske et al. [24,25], MacMillan 
[40], Mynemi [44]); methods of reduction to a nonlinear integral equation (Alobaidi [1], 
Kwok [38], Mallier et al. [41,42], Sevcovic [46], Stamicar et al. [50]). In recent papers 
[55, 56], Zhu derived a closed form of the free boundary position in terms of on infinite 
series. We also refer to a recent survey paper by Chadam [12] focusing on free boundary 
problems in mathematical finance. 

We remind ourselves that in the case of a constant volatility there are, in principle, two 
ways how to solve numerically the free boundary problem for the value of an American call 
resp. put option and the position of the early exercise boundary. The first class of algorithms 
is based on reformulation of the problem in terms of a variational inequality (see Kwok [38] 
and references therein). The variational inequality can be then solved numerically by the 
so-called Projected Super Over Relaxation method (PSOR for short). An advantage of this 
method is that it gives us immediately the value of a solution. A disadvantage is that one 
has to solve large systems of linear equations iteratively taking into account the obstacle for 
a solution, and, secondly, the free boundary position should be deduced from the solution a 
posteriori. Moreover, the PSOR method is not directly applicable for solving the problem 
©-(O when the diffusion coefficient a may depend on the second derivative of a solution 
itself. The second class of methods is based on derivation of a nonlinear integral equation 
for the position of the free boundary without the need of knowing the option price itself (see 
e.g. Evans, Kuske-Keller [21,36], Mallier and Alobaidi [1,41,42], Sevcovic et al. [46,50], 
Chadam et al. [12, 13]. In this approach an advantage is that only a single equation for the 
free boundary has to be solved provided that a is constant; a disadvantage is that the method 
is based on integral transformation techniques and therefore the assumption a is constant is 
crucial. 

Higher order finite difference approximations of the free boundary problems for call 
or put options are discussed in recent papers by Ankudinova and Ehrhardt [2, 3], Ehrhardt 
and Mickens [18] and Zhao et al. [54]. Other interesting analytical and numerical methods 
for evaluating the early exercise boundary have been recently studied by e.g. Milstein 
et al. [43] (Monte-Carlo methods), Widdicks et al. [52] (singular pertubation techniques 
and asymptotic expansions), Cho et al. [14] (parameter estimation methods), Grandits and 
Schachinger et al. [26] (tracking of a discontinuity method), Imai et al. [31] (numerical 
method for generalized Lelend model). 

In this survey we recall an iterative numerical algorithm for solving the free boundary 
problem for an American type of options in the case the volatility a may depend on the 
option and asset values as well as on the time T — t to expiry as well as for American type 
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of Asian options. A key idea of this method consists in transformation of the free boundary 
problem into a semilinear parabolic equation defined on a fixed spatial domain coupled 
with a nonlocal algebraic constraint equation for the free boundary position. It has been 
proposed and analyzes by the author in a series of papers [46-48,50]. Since the resulting 
parabolic equation contains a strong convective term we make use of the operator splitting 
method in order to overcome numerical difficulties. A full space-time discretization of 
the problem leads to a system of semi-linear algebraic equations that can be solved by an 
iterative procedure at each time level. 

The rest of the chapter has the following organization: in the next section we recall the 
well known nonlinear generalization of Black-Scholes equations due to Frey and Stremme, 
Barles and Soner and Kratka, Jandacka and Sevcovic, resp. We also present qualitative 
and quantitative properties of the nonlinear models for pricing European style of options 
with special focus on the Risk adjusted pricing methodology (RAPM) due to Kratka [37], 
Jandacka and Sevcovic [32]. Section 3 is devoted to the free boundary problem for pricing 
American options by means of a linear Black-Scholes equation, i.e. a > is constant. 
This is important in order to understand important steps of the fixed domain transformation 
method. A resulting system of transformed equations consists of a nonlocal parabolic equa- 
tion defined on a fixed domain with time depending coefficients and an algebraic constraint 
equation for the free boundary position. Since the volatility a is constant in this case, by 
using Sine and Cosine Fourier transformations the system of transformed equations can be 
further simplified and reduced to a single nonlinear integral equation for the free boundary 
position. We show how this integral equation can be utilized in order to obtain qualitative 
properties of the free boundary position (early exercise behavior, long time behavior) as 
well as quantitative properties (a fast and stable numerical scheme for computing the free 
boundary function). In Section 4 we discuss a transformation method applied to a class of 
nonlinear Black-Scholes equations. We are able to derive a similar system of transformed 
parabolic-algebraic equation to the one from Section 3. However, as the volatility is no 
longer a constant and it may depend on the solution itself the resulting system of equa- 
tions can not be reduced to a single integral equation for the free boundary and it has to 
be solved numerically. We propose a numerical method based on the finite difference ap- 
proximation combined with an operator splitting technique for numerical approximation of 
the solution and computation of the free boundary position. Several numerical results for 
nonlinear Black-Scholes equations with volatility functions a defined as in ((6]) and ([8]) are 
presented. We also compare our methodology with well-known methods for evaluation of 
approximation to the free boundary position in the case the volatility a is constant. We an- 
alyze dependence of the free boundary position with respect to various parameters entering 
expressions © and (HJ. Finally, Section 5 is devoted to a recent application of the trans- 
formation method in the case of American style of Asian options in which the strike price 
is an arithmetical average of underlying asset prices. Although the volatility a is assumed 
to be constant, due to a specific character of the Black-Scholes equation for pricing Asian 
option the resulting transformed system of equation cannot be further reduced and has to be 
solved numerically by a slight modification of the numerical method discussed in Section 
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4. We finish the last section by presentation of several illustrative examples describing the 
early exercise boundary for American style of Asian call options with floating strike. 

2 Risk adjusted methodology model 

The aim of this section is to present one of nonlinear generalizations of the classical Black- 
Scholes equation with a volatility a of the form (O in a more detail. We focus on the 
so-called Risk adjusted pricing methodology model due to Kratka [37] and straightforward 
generalization by Jandacka and Sevcovic [32] (see also [47]). In this model both the risk 
arising from nontrivial transaction costs as well as the risk from unprotected volatile port- 
folio are taken into account. Their sum representing the total risk is subject of minimiza- 
tion. The original model was proposed by [37]. In [32] we modified Kratka's approach 
by considering a different measure for risk arising from unprotected portfolio in order to 
construct a model which is scale invariant and mathematically well posed. These two im- 
portant features were missing in the original model of Kratka. The model is based on the 
Black-Scholes parabolic PDE in which transaction costs are described by the Hoggard, 
Whalley and Wilmott extension of the Leland model (cf. [29,30,38]) whereas the risk from 
a volatile portfolio is described by the average value of the variance of the synthesized 
portfolio. Transaction costs as well as the volatile portfolio risk depend on the time-lag 
between two consecutive transactions. We define the total risk premium as a sum of trans- 
action costs and the risk cost from the unprotected volatile portfolio. By minimizing the 
total risk premium functional we obtain the optimal length of the hedge interval. 

Concerning the dynamics of an underlying asset we will assume that the asset price 
S = S(t),t > 0, follows a geometric Brownian motion (Q]) with a drift p, standard deviation 
a > and it may pay continuous dividends, i.e. dS = (p — q)Sdt + aSdW where dW 
denotes the differential of the standard Wiener process and q > is a continuous dividend 
yield rate. This assumption is usually made when deriving the classical Black-Scholes 
equation (see e.g. [30, 38]). Similarly as in the derivation of the classical Black-Scholes 
equation we construct a synthesized portfolio II consisting of a one option with a price V 
and 5 assets with a price S per one asset: 

U = V + 5S. (11) 

We recall that the key idea in the Black-Scholes theory is to examine the differential All 
of equation (fTTb . The right-hand side of (fTTT) can be differentiated by using Ito's formula 
whereas portfolio's increment AH(t) = Il(t + At) — H(t) of the left-hand side can be 
expressed as follows: 

An = rUAt + 5qSAt (12) 

where r > is a risk-free interest rate of a zero-coupon bond. In the real world, such a 
simplified assumption is not satisfied and a new term measuring the total risk should be 
added to (fl2l) . More precisely, the change of the portfolio II is composed of two parts: 
the risk-free interest rate part rflAt and the total risk premium: r^SAi where is a risk 
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premium per unit asset price. We consider a short positioned call option. Therefore the 
writer of an option is exposed to this total risk. Hence we are going to price the higher Ask 
option price - an offer to sell an option. It means that All = rllAi — rpSAt. The total 
risk premium r p consists of the transaction risk premium rxc and the portfolio volatility 
risk premium ryp, i.e. rp = r?c + fVP- Hence 



An = rllAt + 5qS At -(r TC + r yP ) S At. (13) 

Our next goal is to show how these risk premium measures rxci r VP depend on the time lag 
and other quantities, like e.g. a, S, V, and derivatives of V. The problem can be decomposed 
in two parts: modeling the transaction costs measure r^c an d volatile portfolio risk measure 
rvP- 

We begin with modeling transaction costs. In practice, we have to adjust our portfolio 
by frequent buying and selling of assets. In the presence of nontrivial transaction costs, 
continuous portfolio adjustments may lead to infinite total transaction costs. A natural way 
how to consider transaction costs within the frame of the Black-Scholes theory is to follow 
the well known Leland approach extended by Hoggard, Whalley and Wilmott (cf. [29,38]). 
We will recall an idea how to incorporate the effect of transaction costs into the governing 
equation. More precisely, we will derive the coefficient of transaction costs rye occurring 
in (TT3T ). Let us denote by C the round trip transaction cost per unit dollar of transaction. 
Then 

C = (S ask - S bid )/S (14) 

where S as k and Sbid are the so-called Ask and Bid prices of the asset, i.e. the market price 
offers for selling and buying assets, resp. Here S = (S as f. + Sind)/2 denotes the mid value 
of the underlying asset price. 

In order to derive the term vpc m 03 measuring transaction costs we will assume, 
for a moment, that there is no risk from volatile portfolio, i.e. ryp = 0. Then AV + 
5 AS = An = rllAt + 5qSAt + r TC SAt. Following Leland's approach (cf. [29]), using 
Ito's formula and assuming <5-hedging of a synthetised portfolio n one can derive that the 
coefficient rpc of transaction costs is given by the formula: 

r TC = ^\d 2 s V\^= (15) 

(see [29, Eq. (3)]). It leads to the well known Leland generalization of the Black-Scholes 
equation © in which the diffusion coefficient is given by (@]) (see [29, 39] for details). 

Next we focus our attention to the problem how to incorporate a risk from a volatile 
portfolio into the model. In the case when a portfolio consisting of options and assets is 
highly volatile, an investor usually asks for a price compensation. Notice that exposure to 
risk is higher when the time-lag between portfolio adjustments is higher. We shall propose a 
measure of such a risk based on the volatility of a fluctuating portfolio. It can be measured 
by the variance of relative increments of the replicating portfolio U = V + 5S, i.e. by 
the term var((AU)/ S). Hence it is reasonable to define the measure ryp of the portfolio 
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volatility risk as follows: 

var (4P) 

r V p = R jf- 1 - (16) 

In other words, ryp is proportional to the variance of a relative change of a portfolio per 
time interval At. A constant R represents the so-called risk premium coefficient. It can 
be interpreted as the marginal value of investor's exposure to a risk. If we apply Ito's 
formula to the differential An = AV + 5 AS we obtain An = (d s V + 5) aSAW + 
\a 2 S 2 T{AW) 2 + g where T = d 2 s V and Q = (d s V + 5)pSAt + d t VAt is a deterministic 
term, i.e E(Q) = Q in the lowest order At - term approximation. Thus 

An - £(An) = (d s V + 5) aScpVAt + ^a 2 S 2 {4> 2 - l)TAt 

where 4> is a random variable with the standard normal distribution such that AW = <p\/~Ai. 
Hence the variance of An can be computed as follows: 

war(AII) = £ ([All - E{ATl)} 2 ) = E (\{d s V + 5)aS<j>y/M + \a 2 S 2 T (0 2 - l) At]' 

Similarly, as in the derivation of the transaction costs measure rxc we assume the ^-hedging 
of portfolio adjustments, i.e. we choose 5 = —dsV. Since E((cf) 2 — l) 2 ) = 2 we obtain an 
expression for the risk premium ry p in the form: 

r VP = ^Ra 4 S 2 T 2 At. (17) 

Notice that in our approach the increase in the time-lag At between consecutive transac- 
tions leads to a linear increase of the risk from a volatile portfolio where the coefficient 
of proportionality depends the asset price S, option's Gamma, T = dgV, as well as the 
constant historical volatility a and the risk premium coefficient R. 



2.1 Risk adjusted Black-Scholes equation 

The total risk premium rp = rpc + r VP consists of two parts: transaction costs premium 
rpc an d the risk from a volatile portfolio ryp premium defined as in (031) and (fTTT ). resp. 
We assume that an investor is risk aversive and he/she wants to minimize the value of the 
total risk premium rp. For this purpose one has to choose the optimal time-lag At between 
two consecutive portfolio adjustments. As both rpc as wei l as r VP depend on the time-lag 
At so does the total risk premium rp. In order to find the optimal value of At we have to 
minimize the following function: 

(j\Y\frS 1 1 

At r R = r TC + r V p = 7= + -Ra 4 S 2 T 2 At . 

V2tt ^At 2 

The unique minimum of the function At rp(At) is attained at the time-lag At op t = 
K 2 /(a 2 \Sr\ 3) where K = {C/{Ry/2n)s. For the minimal value of the function At ^ 
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(18) 



Taking into account both transaction costs as well as risk from a volatile portfolio effects 
we have shown that the equation for the change All of a portfolio II read as: 

AV + 5 AS = AUAt = rUAt + 5qSAt - r R SAt 

where r R represents the total risk premium, r R = rxc + f"VP- Applying Ito's lemma 
to a smooth function V = V(S, t) and assuming the <5-hedging strategy for the portfolio 
adjustments we finally obtain the following generalization of the Black-Scholes equation 
for valuing options: 



8V a 2 2 d 2 V 

~m + T ds 2 " + (r 



dV 

q)s ds 



I'V — trS = . 



By taking the optimal value of the total risk coefficient r R derived as in (1181 ) the option price 
V is a solution to the following nonlinear parabolic equation: 
(Risk adjusted Black-Scholes equation) 



d 2 v 

dS 2 



+ (r 



dV 

q)s ds 



rV = 0. 



(19) 



C 2 R 
P — 



and u p with u 



SdlV and p 



1/3 stands for the signed power 



where pL = 3 

function, i.e. u p = |u| p_i u. In the case there are neither transaction costs (C = 0) nor 
the risk from a volatile portfolio (R = 0) we have \i = 0. Then equation ( fT9l reduces 
to the original Black-Scholes linear parabolic equation Q. We note that equation ( fT9l 

"21 

is a backward parabolic PDE if and only if the function (3(H) = ^-(1 + fiHa)H is an 
increasing function in the variable H := ST 
H > 0. 



SdgV. It is clearly satisfied if \i > and 
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As it is usual in the classical Black-Scholes theory for European style of options (cf. 
[30,38]) we consider the change of independent variables: x := In (J;) , x G R , r := 
T — t , r G (0, T) . As equation ( fl9l ) contains the term ST = Sd 2 s V it is convenient to 
introduce the following transformation: H{x,t) := ST = SdgV(S,t). Now we are in 
position to derive an equation for the function H. It turns out that the function H{x,t) is 
a solution to a nonlinear parabolic equation subject to the initial and boundary conditions. 
More precisely, by taking the second derivative of equation dT9l > with respect to x we obtain, 
after some calculations, that H = H(x, r) is a solution to the quasilinear parabolic equation 

r G (0, T), x G R (see [32]). Henceforth, we will refer to (l20l as a T equation. A solution 
.ff to (l20l) is subjected to the initial condition at r = 0: 

H(x,0)=H(x), xGR, (21) 

where H(x) is the Dirac 5 function H(x) = 5(x). For the purpose of numerical approxi- 
mation we approximate the initial Dirac delta function by H(x) = N'(d)/(ay/r*) where 
r* > is sufficiently small, N(d) is the cumulative distribution function of the normal dis- 
tribution, and d = {x + (r — q — a 2 /2)r* / ay/r*). It corresponds to the value H = SdgV 
of a call (put) option valued by a linear Black-Scholes equation with a constant volatility 
a > at the time T — r* close to expiry T when the time parameter < r* <C 1 is 
sufficiently small. In the case of call or put options the function H is subjected to boundary 
conditions at x = ±oo, 

H(-oo,t) =H(oo,t) =0, t€(0,T). (22) 



A numerical discretization scheme of the T equation ([201 based on finite volume approxi- 
mation has been discussed in [32] in more details. 



2.2 Pricing of European style of options by the RAPM model 

Let us denote V(S, t; C, a, R) the value of a solution to (|T9T > with parameters C, a, R. Sup- 
pose that the coefficient of transaction costs C is known from and is given by (THl) . In real 
option market data we can observe different Bid and Ask prices for an option, Vud < Vask-> 
resp. Let us denote by V m id the mid value, i.e. V m id = \(Yud + Vask)- By the RAPM 
model we are able to explain such a Bid- Ask spread in option prices. The higher Ask price 
corresponds to a solution to the RAPM model with some nontrivial risk premium R > 
and C > whereas the mid value V m id corresponds to a solution V(S, t) for vanishing risk 
premium R = 0, i.e. to a solution of the linear Black-Scholes equation ([2]). An illustrative 
example of Bid- Ask spreads captured by the RAPM model is shown in Fig. [2l 

To calibrate the RAPM model we seek for a couple (o-rapm, R) of implied RAPM 
volatility grapm an d risk premium R such that V as k = V(S,t;C, cfrapm > R) and V m id = 
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Figure 2: A comparison of Bid and Ask option prices computed by means of the RAPM model. 
The middle dotted line is the option price computed from the Black-Scholes equation. We chose 
a = 0.3, /i = 0.2, r = 0.011, E = 25 and T = 1 (left) and T = 0.3 (right). 




Figure 3: Intra-day behavior of Microsoft stocks (April 4, 2003) and shortly expiring call options 
with expiry date April 19, 2003. Computed implied volatilities (Jrapm and risk premium coeffi- 
cients R. 

V(S,t;C,aRAPM,0). Such a system of nonlinear equation for cfrapm and R can be 
solved by by means of the Newton-Kantorovich iterative method (cf. [32]). 

As an example we considered sample data sets for call options on Microsoft stocks. We 
considered a flat interest rate r = 0.02, a constant transaction cost coefficient C = 0.01 
estimated from (fT4l . and we assumed that the underlying asset pays no dividends, i.e. q = 0. 
In Fig. [3] we present results of calibration of implied couple {<trapm,R)- Interestingly 
enough, two call options with higher strike prices E = 25, 30 had almost constant implied 
risk premium R. On the other the risk premium of an option with lowest E = 23 was 
fluctuating and it had highest average of R. 

Finally, in Fig. |we present one week behavior of implied volatilities and risk premium 
coefficients for the Microsoft call option on E = 25 expiring at T = April 19, 2003. In 
the beginning of the investigated period the risk premium coefficient R was rather high and 
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Figure 4: One week behavior of Microsoft stocks (March 20 - 27, 2003) and call options with 
expiration date April 19, 2003. Computed implied volatilities <trapm and risk premiums R. 



fluctuating. On the other hand, it tends to a flat value of R 5 at the end of the week. 
Interesting feature can be observed at the end of the second day when both stock and option 
prices went suddenly down. The time series analysis of the implied volatility cfrapm from 
first two days was unable to predict such a behavior. On the other, high fluctuation in the 
implied risk premium R during first two days can send a signal to an investor that sudden 
changes can be expected in the near future. 



3 Transformation method for a linear Black-Scholes equation 



One of the interesting problems in this field is the analysis of the early exercise boundary 
and the optimal stopping time for American options on stocks paying a continuous dividend. 
It can be easily reduced to a problem of solving a certain free boundary problem for the 
Black-Scholes equation (cf. Black & Scholes [7]). However, the exact analytical expression 
for the free boundary profile is not known yet. Many authors have investigated various 
approximate models leading to approximate expressions for valuing American call and put 
options (see e.g. [24, 25, 33, 35, 36, 44, 45] and recent papers by Zhu [55, 56], Alobaidi et 
al. [1,41,42], Stamicar et al. [50] and the survey paper by Chadam [12] and other references 
therein). For the purpose of studying the free boundary profile near expiry, many different 
integral equations have been derived [6, 11,30,40]. 

Let us recall that the equation governing the time evolution of the price V(S, t) of the 
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American call option is the following parabolic PDE: 



V(0, t) = 0, V(S f (t),t) = S f (t) - E , —(S f (t),t) = 1 , 
V(S,T) = max(S - E, 0) , 




< S < S f (t), 0<t<T, 



(23) 



defined on a time-dependent domain S G (0,Sf(t)), where t 6 (0,T). As usual S 1 > 
is the stock price, E > is the exercise price, r > is the risk-free rate, g > is the 
continuous stock dividend rate and 



is a constant volatility of the underlying stock process. 

The main purpose of this section is to present an alternative integral equation which 
will provide an accurate numerical method for calculating the early exercise boundary near 
expiry. The derivation of the nonlinear integral equation is based on the Fourier transform. 
A solution to this integral equation is the free boundary profile. The novelty of this approach 
consists in three steps: 

1 . The fixed domain transformation. 

2. Derivation of a parabolic PDE for the so-called synthetic portfolio with a nonlin- 
ear nonlocal constraint between the free boundary position and the solution of this 
parabolic equation itself. 

3. Construction of a solution by means of Fourier sine and cosine integral transforms. 

Throughout this section we restrict our attention to the case when r > q > 0. It is well 
known that, for r > q > 0, the free boundary q{t) = St(T — r) starts at q(0) = rE/q, 
whereas g(0) = E for the case r < q (cf. Dewynne et al. [16], Kwok [38]). Thus, the 
free boundary profile develops an initial jump in the case r > q > 0. Notice that the case 
< r < q can be also treated by other methods based on integral equations. Kwok [38] 
derived another integral equation which covers both cases < r < q, as well as r > q > 
(see equation (145T ) and Remark I3.2I ). However, in the latter case equation (|43T > becomes 
singular as t — > T~ , leading to numerical instabilities near expiry. 

In the rest of this section to investigate the behavior of the free boundary Sf(t). We 
present a method developed by Sevcovic in [46] of reducing the free boundary problem 
for (l23l to a nonlinear integral equation with a singular kernel. Notice that our method of 
reducing the free boundary problem to a nonlinear integral equation can be also successfully 
adopted for valuing the American put option paying no dividends (cf. Stamicar et al. [50]). 



a = const > 
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3.1 Fixed domain transformation of the free boundary problem 

In this section we will perform a fixed domain transformation of the free boundary problem 
(|23T ) into a parabolic equation defined on a fixed spatial domain. As it will be shown below, 
imposing of the free boundary condition will result in a nonlinear time-dependent term 
involved in the resulting equation. To transform equation d23l > defined on a time dependent 
spatial domain (0, Sf(t)), we introduce the following change of variables: 

T = T-t, jc = ln^^ where q(t) = S f (T-r). (24) 

Clearly, r G (0, T) and x G (0, oo) whenever S G (0, Sf(t)). Let us furthermore define the 
auxiliary function n = II(x, r) as follows: 

U(x,r) = V(S,t)-S^(S,t). (25) 

Notice that the quantity II has an important financial meaning as it is a synthetic portfolio 
consisting of one long option and A = S underlying short positioned stocks. It follows 
from (O that 

dn n2 d 2 V <9 2 n dU n ,d 3 V 
= |_ 2 = —S 

dx dS 2 ' dx 2 dx dS 3 ' 

du | gdu _ s d 2 v av 

dr q dx dSdt dt ' 

where g = dg/dr. Now assuming that V = V(S,t) is a sufficiently smooth solution of 
(|23T ). we may differentiate (T23T ) with respect to S. Plugging expressions (l26l) into (|23T ). we 
finally obtain that the function II = Ii(x, r) is a solution of the parabolic equation 

dU , 811 a 2 d 2 U 

* + °< T) &-Ta? + rn = ' <27) 

x G (0,oo), r G (0, T), with a time-dependent coefficient 

. , q(t) a 2 

It follows from the boundary condition V(S f (t),t) = S f (t) - E and V s (S f (t),t) = 1 that 

11(0, r) = —E, n(+oo, r) = . (28) 

The initial condition II(a;,0) can be deduced from the pay-off diagram for V(S, T). We 
obtain 

n(x,o) = j - E (29 ) 

otherwise. 
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Notice that equation (|27T ) is a parabolic PDE with a time-dependent coefficient a(r). In 
what follows, we will show the function a(r) depends upon a solution II itself. This depen- 
dence is non-local in the spatial variable x. Moreover, the initial position of the interface 
g(0) enters the initial condition H(x, 0). Therefore we have to determine the relationship 
between the solution Ii(x, r) and the free boundary function g(r) first. To this end, we 
make use of the boundary condition imposed on V at the interface S = Sf(t). Since 
S f (t) -E = V(S f (t),t) we have 

d _ . . dV , _ . . . d _ . . <9F , _ , , . 
^) = ^(5 / (t),t)-5 / (t) + -(5 / (t),t). 

As ^(Sf(t),t) = 1 we obtain %{S,t) = at 5 = S/(i). Assuming the function 
n^. has a trace at x = 0, and taking into account d26l ), we may conclude that, for any 

t = T-r€ [0,T), 



.'2 



.,3 2 y an av" 



If ^(5,t) -» %{Sf(t),t) = as 5 -> 5/(t) - , then it follows from the Black-Scholes 
equation d23l that 

(j 2 an 

(r - </)f?(r) + y^(0,r) - r(<?(r) - E) 

I BV f)V rr 2 rP-V \ 

As a consequence we obtain a nonlocal algebraic constraint between the free boundary 
function q{t) and the boundary trace d x ll(0, r) of the derivative of the solution II itself: 

rE a 2 dU, 

£ T = _+ 0,r forO<r<T. (30) 

q 2q ox 

It remains to determine the initial position of the interface £»(0). According to Dewynne et 
at. [16] (see also Kwok [38]), the initial position g(0) of the free boundary is rE/q for the 
case < q < r. Alternatively, we can derive this condition from (|27T>— d29l> by assuming the 
smoothness of II in the x variable up to the boundary x = uniformly for r — > + . More 
precisely, we assume that 

an , . , an. , , an . . 

lim — — (0, r)= lim — — (x,t) = lim — — (x,0)=0 
r^0+ OX t~>0+,x~>0+ OX x->0+ ox 

because H(x, 0) = — E for x close to + . By (l30l) we obtain 

rE 

Q(0) = — • (3D 
Q 
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In summary, we have shown that, under suitable regularity assumptions imposed on a solu- 
tion n to (|27T ). (T28T ). (l29l) . the free boundary problem (1231) can be transformed into the initial 
boundary value problem for parabolic PDE 



dU a(r) dU a 2 d 2 n (| _ () 
dr dx 2 Ore 2 

n(0, r) = n(+oo, r) = 0, x > 0, r G (0, T), (32) 



n(x,o) 



— E for x < In (r/q) 
otherwise, 



where a(r) = 4j4 + r — q — ^ and 



a 2 an, N , , r.E 
e (r) = — + — — 0,r), e(0) = — . (33) 
q 2q ox q 



We emphasize that the problem (1321 ) constitutes a nonlinear parabolic equation with a non- 
local constraint given by (l33l . 



Remark 3.1 /« our derivation of the free boundary function q(t) and its initial condition 
q(0) we did not assume that the solution V(S, t) is C 2 smooth up to the free boundary S = 
Sf(t). Such an assumption would lead to an obvious contradiction F := d 2 V/dS 2 = at 
S = Sf(t). On the other hand, the jump in F at the free boundary is the only driving force 
for the evolution of the free boundary function q (see (I30D ). Construction of a PDE for the 
synthetic portfolio function II is a crucial step in our approach because the derivative d x FL 
admits a trace at the boundary x = and the unknown free boundary function q can be 
determined via (13 3i 



3.2 Reduction to a nonlinear integral equation 

The main purpose of this section is to show how the fully nonlinear nonlocal problem (|32~I)- 
(l33l can reduced to a single nonlinear integral equation for q(t) giving the explicit formula 
for the solution II(x, r) to (|32~1) . The idea is to apply the Fourier sine and cosine integral 
transforms (cf. Stein & Weiss [51]). Let us recall that for any Lebesgue integrable function 
/ G the sine and cosine transformations are defined as follows: J r s{f){(jj) = 

fo° f( x ) sinwx dx, J r c(/)(w) = / °° f(x) cos ujx dx. Their inverse transforms are given 
by J r g 1 (g)(x) = ^ / °° g(u) sinuixdu, J : ^, l (g)(x) = \ g(u) cosuxduj. Now we 
suppose that the function q(t) and subsequently a(r) are already know. Let II = H(x, r) 
be a solution of (l32l ) corresponding to a given function a(r). Let us denote 

p(w,r)=^ s (n(.,r))(a;) > q(u, r) = ^ c (n(., T ))(u) (34) 

where lo G M + ,r G (0, T). By applying the sine and cosine integral transforms, taking 
into account their basic properties and (l33l) . we finally obtain a linear non-autonomous lo- 
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parameterized system of ODEs 



4^p{u, T ) ~ a(r)ujq(uj, r) + a(w)p(w, r) = -Buy, 

-^-q(u), t) + a(r)up(u, r) + a(u)q(u, r) = —EaOr) - qg(j) + rE, (35) 
ctr 

for the sine and cosine transforms of II (see d34l) ) where 

a(w) = \{o- 2 oo 2 + 2r). 

The system of equations (|33T > is subject to initial conditions at r = 0, p(o>,0) = 
Ts(H{-, 0)(w), <7(u;, 0) = J r c(n(., 0)(w). In the case of the call option, we deduce from 
the initial condition for II (see d32l >) that 



E 



p(uj, 0) = — I cos I uln - 1 — 1 



E ( r 

q(oJ, 0) = sin a; In - 

uj \ q 



(36) 



Taking into account d36l > and by using the variation of constants formula for solving linear 
non-autonomous ODEs, we obtain an explicit formula for p(lu,t) = —Elo^ 1 + p(lu,t) 
where 

E -a(u)T 



p( W) r) = -e- a ^ T cos(w(A(r,0) + \n(r/q))) 



-a(u>)(T—s) 



rE 



cos(uM(r, s)) + (rE — qg(s)) sin(uM(T, s)) 



ds. (37) 



Here we have denoted by A the function defined as 



A(t,s) 



Q{r) 



a- 



(38) 



As F s l {u) l ) = 1 we have II(a;, r) = J t s 1 (p(oj,t)) = —E + ^ p(u, r) sin(a;x) duj. 
From (l33l we conclude that the free boundary function g satisfies the following equation: 



rE a 2 f°° 
g(r) = 1 / ujp(ui , t) du 

q Jo 



(39) 



Taking into account (I37T ) and (l39l we end up with the following nonlinear singular integral 
equation for the free boundary function g(r): 



g(r) ~ rE (l I ° expf rr + Mr/g)) s 

8[ ] ~ q{ + r^hTr P V 2<rV 



+ 



1 



2tt Jo 



cr \ rE 1 / r — s 



exp f _ r ( r _ s ) _ _^Zi£l_ 



T — S 



(40) 



(is , 
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where the function A depends upon g via equation (138T ). To simplify this integral equation, 
we introduce a new auxiliary function H : [0, vT] — > R as follows: 

g(r) = — (l + aV2H(V^)) . (41) 



Using the change of variables s = £ 2 cos 2 9, one can rewrite the integral equation d39l) in 
terms of the function if as follows: 

7T 

# (0 =/ff(0 + 4= / 2 K cos ^ - 2cotg £T(£ cos % H (£, *)] e"^ sin2 d0, 
V 71 " Jo 

(42) 

where 

* -JL-lnf 1 + ) + ^ sm Q (43) 

for£e [0,VT], 0e (0,7r/2), 

cr 2 

and 



Notice that equations (1401 ) and (1421 are integral equations with a singular kernel (cf . Gripen- 
berg et al. [27]). 

Remark 3.2 Kwok [38] derived another integral equation for the early exercise boundary 
for the American call option on a stock paying continuous dividend. According to Kwok [ 38, 
Section 4.2.3], q(t) satisfies the integral equation 

q(t) =E + e {r)e~ qT N(d) - Ee~ rr N{d - ay/r) 

+ f qg{T)e- q ^N(d ( ) - rEe~ ri N(d ( - o^fl)d£ (45) 
Jo 

where 

a^F \ E J ' 5 cV£ \q(t-0J 

and N(u) is the cumulative distribution function for the normal distribution. The above 
integral equation covers both cases: r < q as well as r > q. However, in the case r > q 
this equation becomes singular as r — > + . 

In the rest of this section we derive a formula for pricing American call options based 
on the solution q to the integral equation (|42)) . With regard to (l25l) . we have 



^ (s- l v{s,t)) = -s- 2 u (in (s-yr -t)),T-t) 
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Taking into account the boundary condition V(Sf(t),t) = Sf(t) — E and integrating the 
above equation from S to Sf(t) = g(T — t), we obtain 



V(S,T-t) 



s 

W) 



q(t) -E + 



In 



e(r) 



e x U(x, r)dx) . 



(46) 



It is straightforward to verify that V given by (1461 is indeed a solution to the free bound- 
ary problem (|23T ). Inserting the expressions (l37l) and (l39l) (recall that H(x, r) = — E + 
§ Jo° Pi^i T ) smwa:; ^-0 i nto d46l) . we end up with the formula for pricing the American 
call option: 



V(S,T-r) 



+ 



S-E+ -^-E I 2 (A(r, 0) + Hr/q)MQ(r)/S), r) 
S 



Q{r) 



rE I 2 (A(t,8)Mq(t)/S),t - s) 



(47) 



+{rE - qg(s)) h(A(r, s)Mq(t)/S),t - s) 



for any S G (0, S f (t)) and t G [0, T], where A(r, s) = In jfej + (r - qf- £)(r - s). Here 



-(r-<x 2 /2)r 



-A - a 2 r L 



e~ rT e L 



M 



A-L 2L 



e~ A M 



A — a 2 r L 



aV2r ' aV2r 



-(r-a 2 /2)r 



e A M 



A - cj 2 r L 



(48) 



and 



M(x, y) = erf (a; + y) — erf(x) 



+ e~ A M 



x+y 



A — a t L 



We will refer to (l47l as the semi-explicit formula for pricing the American call option. 
We use the term 'semi-explicit' because (l47l) contains the free boundary function q{t) = 
Sf(T — r) which has to be determined first by solving the nonlinear integral equation (|42l . 



3.3 Numerical experiments 

In this section we focus on numerical experiments. We will compute the free boundary 
profile 

rE 



S f (t) = g(T-t) 



1 + aV2H(^T - t) 



(49) 



(see d36l >) by solving the nonlinear integral equation (l42l . We also present a comparison of 
the results obtained by our methods to those obtained by other known methods for solving 
the American call option problem. 
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The computation of a solution of the nonlinear integral equation is based on an iterative 
method. We will construct a sequence of approximate solutions to ((42l . Let H be an initial 
approximation of a solution to (l42l . If we suppose H°(£) = h\£, then, plugging this ansatz 
into (|42l yields the well-known first order approximation of a solution i/(£) in the form 

F°(£) = 0.451381 £ 

i.e. q°(t) = — (1 + 0.638349 (Jy/r) (cf. Sevcovic [46]). This asymptotics is in agreement 
with that of Dewynne et al. [16]. For n = 0, 1, ... we will define the n + 1 approximation 
/yn+i as f n ows: 

i? n+l (e) = 

7T 

/#«(£) + -L / [^cos0-2cotg^ J ff n (ecos0)(7 // n( ( e,e)]e- rf2sin29 ~^«. e )^ (50) 
V 71 " Jo 

for £ G [0, \/r]. With regard to (@3) and (|44]l . we have 
9H«{C,0) = — t=— — —In 



cos 




-re a -(fffln(e,f)+|^i n (i)) J 



2r0F£ 

Notice that the function gu is bounded provided that H is nonnegative and Lipschitz con- 
tinuous on [0, VT}. Recall that we have assumed r > q > 0. Then the function £ /#(£) 
is bounded for £ G [0, VT] and it vanishes at £ = 0. Moreover, if H is smooth then fjj 
is a flat function at £ = 0, i.e. /#(£) = o(£ n ) as £ — > + for all n G N. From the 
numerical point of view such a flat function can be omitted from computations. For small 
values of 6 we approximate the function cotg 9 gH n (£, Q) by its limit as 9 — > + . It yields 
the approximation of the singular term in (l50l > in the form 



co,g W )~l i + fJ^ ) + A { te0 <»«l, 

where i7 = i/ n . In the following numerical simulations we chose e 10 -5 . 

In what follows, we present several computational examples. In Fig. [5] we show five iter- 
ates of the free boundary function Sf(t), where the auxiliary function if(£) is constructed 
by means of successive iterations of the nonlinear integral operator defined by the right- 
hand side of (l42l . This sequence converges to a fixed point of such a map, i.e. to a solution 
of (|42~1) . Parameter values were chosen as E = 10, r = 0.1, ex = 0.2, q = 0.05, T = 1. 
Fig. |5]depicts the final tenth iteration of approximation of the function Sf(t). The mesh 
contained 100 grid points. One can observe very rapid convergence of iterates to a fixed 
point. In practice, no more than six iterates are sufficient to obtain the fixed point of ((42l . 
It is worth noting that in all our numerical simulations the convergence was monotone, i.e. 
the curve moves only up in the iteration process. 
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Figure 6: Long-time behavior of Sf(t). 



In Fig. [6] we show the long time behavior of the free boundary Sf(t),t S [0, T], for 
large values of the expiration time T. For the parameter values T = 50, r = 0.1, q = 
0.05, a = 0.35 and E = 10 the theoretical value of 5/(+oo) is 36.8179 (see Dewynne et 
al. [16]). 

In Table[T] we show a comparison of results obtained by our method based on the semi- 
explicit formula (1471 and those obtained by known methods based on trinomial trees (both 
with the depth of the tree equal to 100), finite difference approximation (with 200 spatial 
and time grids) and analytic approximation of Barone-Adesi & Whaley (cf. [6], [30, Ch. 
15, p. 384]), resp. It also turned out that the method based on solving the integral equation 
(l42l is 5-10 times faster than other methods based on trees or finite differences. The reason 
is that the computation of V(S, t) for a wide range of values of S based on the semi-explicit 
pricing formula (1471 is very fast provided that the free boundary function g has already been 
computed. 

In Fig. [7] the early exercise boundary <S/(t) is computed for various values of the pa- 
rameter q. In these computations we chose E = 10, T = 0.01, a = 0.45. Of interest is the 
case where q is close to r (q = 0.09 and r = 0.1). 



Transformation methods for linear and nonlinear Black-Scholes equations 



23 



Table 1 : Comparison of the method based on formula d48b with other methods for the parameter 
values E = 10,T = l,a = 0.2, r = 0.1, q = 0.05. The position 5/(0) = g(T) of the free 
boundary at t = (i.e. at r = T) was computed as S/(0) = g(T) = 22.3754 



Method \ The asset value S 


15 


18 


20 


21 


22.3754 


Our method 


5.15 


8.09 


10.03 


11.01 


12.37 


Trinomial tree 


5.15 


8.09 


10.03 


11.01 


12.37 


Finite differences 


5.49 


8.48 


10.48 


11.48 


12.48 


Analytic approximation 


5.23 


8.10 


10.04 


11.02 


12.38 




Figure 7: The early exercise boundary Sf(t) for r = 0.1, q = 0.05. The early exercise boundary 
S f {t) for r = 0.1, 9 = 0.09. 
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3.4 Early exercise boundary for an American put option 

In this section we present results of transformation method applied to valuation of the early 
exercise boundary for American style of a put option. Recall that the early exercise bound- 
ary problem for American put option can be formulated as follows: 



dV dV a 2 2 d 2 V n^ + ^Tom^o^ 

~dt ~dS ~2 ~dS 2 ~ <t <T, Sf(t) < S < oo , 

dV 

V(+oo, t) = 0, V(S f (t),t) = E- S f (t) , —(S f (t),t) = -1 , (51) 



V(S, T) = max(E - S, 0) , 



defined on a time-dependent domain S G (Sf(t), oo), where t G (0, T) (cf. Kwok [38]). 
Again S > stands for the stock price, E > is the exercise price, r > is the risk-free 
rate and a > is the volatility of the underlying stock process. We shall assume that the 
asset pay no dividends, i.e. q = 0. In order to perform a fixed domain transformation of the 
free boundary problem (1521 we introduce the following change of variables 



S_ 



lnf-^J, where T = T-t,g(r) = S f (T - r) 



Similarly as in the case of a call option we define a synthetised portfolio II for the put 
option U(x,t) = V — SQ^T. Then it is easy to verify that II is a solution to the following 
the parabolic equation 

dU . N dII a 2 d 2 n „ . _ 

U(0,t) = E, n(oo,r)=0, (52) 

— (0,T) = -rE, forrG(0,T), (53) 
n(x, 0) = for x > 0, 

where o(r) = ||^y+r— (see Stamicar a/. [50]). Now, following the same methodology 
of applying the Fourier transform we are able to derive an integral equation for the free 
boundary position. It is just the equation %-^g(0, r) = — rE 1 in which the left hand side is 
expressed by the inverse Fourier transform of the solution II as a weakly singular integral 
depending on the free boundary position g (see [50] for details). We omit the technical 
details here. We just recall that the integral equation for the free boundary function q(t) 
yields the following expression for q(t) 



6 [t) = Ee- {r - 2 T }T e aV2 ^ viT) 



in terms of a new auxiliary function r/(r) (see Stamicar et al. [50] for details). Further 
asymptotic analysis of the integral equation enables us to derive an asymptotic formula for 
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Figure 8: Left: asymptotic approximation vs. binomial method for a = 0.25, r = 0.1, E = 10 
and T — t = 8.76 hrs. (0.001 of a year), MBW approximation vs. the asymptotic solution (T54l > for 
T-t= 0.876 hrs (right). 

t](t) as r — > 0. Namely, 



Next we examine how accurately our asymptotic approximation matches the data from 
the binomial method (cf. Kwok [38]). Near expiry at about one hour, the asymptotic ap- 
proximation matches the data from the binomial method (see Fig. [8]). With a = 0.25, r = 
0.1, E = 10 at 8.76 hours the approximation gives an overestimate but of only 0.4 cents 
(see [50]). We also compared our asymptotic solution with MacMillan, Barone-Adesi, and 
Whaley's [6], [4, 384-386], [40] numerical approximation of the American put free bound- 
ary (see Fig. [8]). They apply a transformation that results in a Cauchy-Euler equation that 
can be solved analytically. For times very close to expiry, one can see that our approxima- 
tion of the free boundary matches the data from the binomial and trinomial methods more 
accurately. 

The numerical and analytical results obtained by the transformation method for solv- 
ing free boundary problem for American put options are in agreement with those obtained 
recently by Zhu [55,56], Mallier et al. [41,42], Kuske and Keller [36], Knessl [34]. 

4 Transformation method for a nonlinear Black-Scholes equa- 
tion 

The main goal of this section is to perform a fixed domain transformation of the free bound- 
ary problem for the nonlinear Black-Scholes equation (|2]) into a parabolic equation defined 
on a fixed spatial domain. For the sake of simplicity we will present a detailed derivation of 
an equation only for the case of an American call option. Derivation of the corresponding 
equation for the American put option is similar. Throughout this section we shall assume 
the volatility a 2 appearing in the Black-Scholes equation to be a function of the option 




as r 



0+. 



(54) 
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price S, time to expiry T — t and S 2 dgV, i.e. 

a = a(S 2 d 2 s V,S,T-t). 
Following (1241) we shall consider the following change of variables: 

T = T-t, x = In (q(t)/S) where q{t) = S f {T - r). 

Then r G (0, T) and x G (0, oo) iff S G (0, Sf(t)). The boundary value x = corresponds 
to the free boundary position 5 = S/(i) whereas x +oo corresponds to the default value 
S = of the underlying asset. Following Stamicar et al. [50] and Sevcovic [46, 47] we 
construct the so-called synthetic portfolio function n = n(x, r) defined as follows: 

U(x,r) = V(S,t)-S^(S,t). (55) 

Again, it represents a synthetic portfolio consisting of one long positioned option and A = 
dsV underlying short stocks. Similarly as in Section 3 we have 

dn ~d 2 V dU gdU d / BV 
— l — I y cj 

dx OS 2 ' 8t q dx dt \ 8S 

where we have denoted g = dg/dr. Assuming sufficient smoothness of a solution 
V = V(S, t) to (O we can deduce from (O a parabolic equation for the synthetic port- 
folio function n = r) 

du „, , i 2 an i d ( 2 an 



ar +« T )-r)fc-5fa^fc / )+ riI - 

defined on a fixed domain x G R, t G (0, T), with a time-dependent coefficient 

&(T) = 44+ r -9 (56) 

and a diffusion coefficient given by: a 2 = a 2 (d x Tl(x, r), g{r)e~ x , r) depending on 
r, x and the gradient d x Ii of a solution n. Now the boundary conditions V(0, t) = 
0,V(S f (t),t) = S f (t) - E and d s V(S f (t),t) = 1 imply 

U(0,t) = -E, n(+oo,r) = 0, 0<r<T, (57) 

and, from the terminal pay-off diagram for V(S, T), we deduce 

n (l , ) = / - E fo "<'» W ,58) 

[ otherwise. 

In order to close up the system of equations that determines the value of a synthetic 
portfolio n we have to construct an equation for the free boundary position q(t). In- 
deed, both the coefficient b as well as the initial condition n(x, 0) depend on the func- 
tion q(t). Similarly as in the case of a constant volatility a (see [46, 50]) we pro- 
ceed as follows: since Sf(t) — E = V(Sf(t),t) and dsV(Sf(t),t) = 1 we have 
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f t S f {t) = d s V(S f (t),t)f t S f (t) + d t V(S f (t),t) and so d t V(S,t) = along the free 
boundary S = Sf(t). Moreover, assuming d x Tl is continuous up to the boundary 
x = we obtain S 2 d 2 s V(S,t) -► d x U(0,r) and Sd s V(S,t) -> e (r) as 5 -> 
Now, by taking the limit S — ► Sf(t)~ in the Black-Scholes equation ((JJ) we obtain 
(r - g)g(r) + ^a 2 d x U(0, r) - r(g(r) - -E) = 0. Therefore 

e(r = — + — a 2 ajio,r ,r — 0,r 

q 2g ox 

for < r < T. The value of g(0) can be easily derived from the smoothness assumption 
made on d x Il at the origin x = 0, r = under the structural assumption 

< q < r (59) 

made on the interest and dividend yield rates r, q (cf. [46,47]). Indeed, continuity of d x Il 
at the origin (0, 0) implies lim T ^ + d x ll(0, r) = d x Tl(0, 0) = lim x _+ + d x Tl(x, 0) = 
because H(x, 0) = —E for x close to + provided \n(r/q) > 0. From the above equation 
for q(t) we deduce g(0) = ^ by taking the limit r — ► + . Putting all the above equations 
together we end up with a closed system of equations for II = Il(x, r) and g = g(r) 

an . . , a 2 , an 1 a , 2 an, 

a7 + Wr)- T )^-^(^) + rn = o, 

n(0, t) = -E, n(+oo, r) = , x > , r G (0, T) , 

n(x,0) = (-f forx<ln(r/g) 
1 ' ' \ otherwise, v ' 

where cr = <T(a x n(x, r), g(r)e~ x , r) , 6(r) = |£4 + r — q and the free boundary position 
q(t) = Sf(T — r) satisfies an implicit algebraic equation 

rE a 2 (d x U(0,T),g(r),T) dU (n rE 
q(t) = — + ^(°' r )' with e (0) = — , (61) 

where r G (0, T). Notice that, in order to guarantee parabolicity of equation (l60l > we have 
to assume that the function p a 2 (p, g{r)e~ x ,r)p is strictly increasing. More precisely, 
we shall assume that there exists a positive constant 7 > such that 

v 2 {p,Z,T)+pd p a 2 (p,Z,T) > 7 >0 (62) 

for any £ > 0, r G (0, T) and pel. Notice that condition ((62]> is satisfied for the RAPM 
model in which a 2 = a 2 (l + fip^£,~^) for any \x > and p > 0. Clearly p = 5 2 aJ.V > 
for the case of plain vanilla call or put options. As far as the Barles and Soner model is 
concerned, we have (1 + ty(a 2 e VT p)) and condition (l62l is again satisfied because 

the function ^ is a positive and increasing function in the Barles and Soner model. 
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Remark 4.1 Following exactly the same argument as in f !46D one can derive an explicit 
expression for the option price V(S, t): 

S ( /- lni ir \ 

V(S, T - t) = -^-r [ q{t) -E + J e x U{x, r)dxj . (63) 



4.1 An iterative algorithm for approximation of the early exercise boundary 

The idea of the iterative numerical algorithm for solving the problem d60l >, d6"T| ) is rather 
simple: we use the backward Euler method of finite differences in order to discretize the 
parabolic equation (l60l > in time. In each time level we find a new approximation of a solu- 
tion pair (IT, q). First we determine a new position of q from the algebraic equation (loTT ). 
We remind ourselves that (even in the case a is constant) the free boundary function q(t) 
behaves like rE/q + 0{t 1 / 2 ) for r -> 0+ (see e.g. [16] or [46]) and so 6(r) = 0(t -1 / 2 ). 
Hence the convective term b(T)d x Il becomes a dominant part of equation (l60l > for small 
values of r. In order to overcome this difficulty we employ the operator splitting technique 
for successive solving of the convective and diffusion parts of equation (l60l) . Since the dif- 
fusion coefficient a 2 depends on the derivative d x Il of a solution II itself we make several 
micro-iterates to find a solution of a system of nonlinear algebraic equations. 

Now we present our algorithm in more details. We restrict the spatial domain x G 
(0, do) to a finite interval of values x G (0, L) where L > is sufficiently large. For prac- 
tical purposes one can take L 3 as it corresponds to the interval S G (Sf(t)e~ L , Sf(t)) 
in the original asset price variable S. The value Sf(t)e~ L is then could be a good approx- 
imation for the default value S = if L 3. Let us denote by k > the time step, 
k = T/m, and, by h > the spatial step, h = L/n where m,n G N stand for the 
number of time and space discretization steps, resp. We denote by Uj an approximation of 
H(xi,Tj), gi « q{tj), V ~ b(rj) where Xj = ih,Tj = jk. We approximate the value of 
the volatility a at the node (xi,Tj) by finite difference as follows: 

4 = al^,W) = a((Iii +1 -lt>/h^^,r j ). 

Then for the Euler backward in time finite difference approximation of equation (l60l) we 
have 

+ (V - -{am d x W - -d x {(a j ) 2 d x lP) + rW = (64) 

and the solution IP = IP (a;) is subject to Dirichlet boundary conditions at x = and 
x = L. We set II (x) = H(x, 0). Now we decompose the above problem into two parts - a 
convection part and a diffusive part by introducing an auxiliary intermediate step IP - 2: 

(Convective part) 

W • - IP'" 1 1 
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{Diffusive part) 

1 ~^-d x W - -d x {{aifd x W) + rW = . (66) 

The idea of the operator splitting technique consists in comparison the sum of solutions to 
convective and diffusion part to a solution of (l64l ). Indeed, if d x TP « d x W~2 then it is 
reasonable to assume that IP computed from the system (I65l)-(l66l) is a good approximation 
of the system (l64l ). 

The convective part can be approximated by an explicit solution to the transport equa- 
tion: 

d T n + b{r)d x fl = for x > 0, r G fa-i, 73] (67) 

subject to the boundary condition 11(0, r) = — E and initial condition Ti(x,Tj-\) = 
IF (a). For American style of call option the free boundary q(t) = Sf(T — t) must 
be an increasing function in r and we have assumed < q < r we have b(r) = 
q(t) I q(t) + r — q > and so prescribing the in-flowing boundary condition 11(0, r) = — E 
is consistent with the transport equation. Let us denote by B(t) the primitive function to 
6(r), i.e. B(t) = In q(t) + (r — q)r. Equation (1671 can be integrated to obtain its explicit 
solution: 

n fa . t) = fn^ x (x - B(t) + Bfa-i)) if x - B(r) + B(Tj-i) > , 
1 ' J \-£ otherwise. 

Thus the spatial approximation IT. 2 can be constructed from the formula 



Tp 2 _ j nJ i f £i = ^ - ln £j +ln^-i - (r - q)k > 0, 

' 1 -E otherwise, 



where a linear approximation between discrete values II? - , i = 0, 1, n, is being used to 
compute the value W~ 1 (x i — In qa + In qj_i — (r — q)k). 

The diffusive part can be solved numerically by means of finite differences. Using 
central finite difference for approximation of the derivative d x lP we obtain 



1 



TP W 5 (VrM 2 TP _TT J 

fc * 2 2/l 

Hence, the vector of discrete values IP = {II?, i = 1, 2, ...,n} at the time level j G 
{1,2, m} satisfies the tridiagonal system of equations 

, • ^/ 1 1 / • "/"If.; "< •■ (70) 
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for i = 1, 2, n, where 



~J = - 'j^ II/; —(g j ) 2 _ A^) 2 (7i) 

ft = pi(e>,W) = l + rk-(ai + i{). 

The initial and boundary conditions at r = and x = 0, L, resp., can be approximated as 
follows: 



n 







-E for Xi < In (r/q) 
for Xi > In (r/q) 



for i = 0, 1, n, and II* = Il£ = 0. 

Next we proceed by approximation of equation (loTT) which introduces a nonlinear con- 
straint condition between the early exercise boundary function q(t) and the trace of the 
solution IT at the boundary x = (S = Sf(t) in the original variable). Taking a finite 
difference approximation of d x TI at the origin x = we obtain 

e> = ?f + Y q a2 ( (n ' " n ° )//j ' ^) ' (72) 



Now, equations d69| ), (1701 and (1721) can be written in an abstract form as a system of 
nonlinear equations: 

gP = F(IP,eP), 
jp-i = T(lPj), (73) 

A(w,^)u j =rp'-5, 

where .F(IP, g 7 ') is the right-hand side of the algebraic equation (1721 , T(IP, ^ ) is the 
transport equation solver given by the right-hand side of d69l ) and A = A(IP , f>?) is a tridi- 
agonal matrix with coefficients given by dTTI ). The system (|73T > can be approximately solved 
by means of successive iterates procedure. We define, for j > 1, IF' = IP , gi' = gp . 
Then the (p + l)-th approximation of IP and gP is obtained as a solution to the system: 

n j-f, P +i = T(IF>^> +1 ), (74) 

.4(iF' p , ^ +1 )rp' p+1 = n J '-5'P +1 . 

Notice that the last equation is a linear tridiagonal equation for the vector IP> P+1 whereas 
gi,p+i- an d jp- 2-P+ 1 can be directly computed from (f72l and (|69l ), resp. If the sequence of 
approximate solutions {(IP ,P , ^' p )}2? =1 converges to some limiting value (IP' 00 , gi' 00 ) as 
p — ► oo then this limit is a solution to a nonlinear system of equations d73l at the time level 
j and we can proceed by computing the approximate solution the next time level j + 1. 
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4.2 Numerical approximations of the early exercise boundary 

In this section we focus on numerical experiments based on the iterative scheme described 
in the previous section. The main purpose is to compute the free boundary profile Sf(t) = 
g(T — t) for different (non)linear Black-Scholes models and for various model parameters. 
A solution (IT, g) has been computed by our iterative algorithm for the following basic 
model parameters: E = 10, T = 1 (one year), r = 0.1 (10% p. a) , q = 0.05 (5% p.a.) and 
a = 0.2. We used n = 750 spatial points and m = 225000 time discretization steps. Such 
a time step k = T/m corresponds to 140 seconds between consecutive time levels when 
expressed in real time scale. In average we needed p < 6 micro-iterates (1741 in order to 
solve the nonlinear system d73l with the precision 10~ 7 . 

4.2.1 Case of a constant volatility - comparison study 

In our first numerical experiment we make attempt to compare our iterative approxima- 
tion scheme for solving the free boundary problem for an American call option to known 
schemes in the case when the volatility a > is constant. We compare our solution 
to the one computed by means of a solution to a nonlinear integral equation for q(t) 
(see also [46, 50]). This comparison can be also considered as a benchmark or test ex- 
ample for which we know a solution that can be computed by a another justified algo- 
rithm. In Fig. |9j part a), we show the function g computed by our iterative algorithm for 
E = 10, T = l,r = 0.1, q = 0.05, a = 0.2. At the expiry T = 1 the value of g(T) was 
computed as: g(T) = 22.321. The corresponding value g(T) computed from the integral 
equation (@0]> (cf. [46]) was g(T) = 22.375. The relative error is less than 0.25%. In the 
part b) we present 7 approximations of the free boundary function g(r) computed for dif- 
ferent mesh sizes h (see Tab. |2] for details). The sequence of approximate free boundaries 
gh,h = hx,h,2,.--, converges monotonically from below to the free boundary function g as 
h I 0. The next part c) of Fig. [9] depicts various solution profiles of a function IL(x, r). 
In order to achieve a reasonable approximation to equation (1721 we need very accurate ap- 
proximation of n(x, r) for x close to the origin 0. The parts d) and e) of Fig. |9]depict the 
contour and 3D plots of the function U(x, r). 

In Tab.|2]we present the numerical error analysis for the distance || g^ — g\\ p measured in 
two different norms (L°° and L 2 ) of a computed free boundary position g^ corresponding 
to the mesh size h and the solution g computed from the integral equation described in (l40l 
(cf. [46]). The time step k has been adjusted to the spatial mesh size h in order to satisfy 
CFL condition a 2 k/h 2 w 1/2. We also computed the experimental order of convergence 
eoc(L p ) for p = 2, oo. Recall that the experimental order of convergence can be defined as 
the ratio: 

eoc(L p ) = ln dl gfc ' ~ gllp) - ln (llg^-i ~ Q\\p) 
In hi — In hi— i 

It can be inteipreted as an exponent a = eoc(L p ) for which we have \\gh — g\\ p = 0(h a ). 
It turns out from Tab.|2]that the conjecture on the order of convergence \\g^ — gW^ = 0(h) 
whereas \\gh — g\\2 = 0(h 3 ^ 2 ) as h — > + could be reasonable. 
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Figure 9: a) A comparison of the free boundary function q(t) computed by the iterative algorithm 
(green solid curve) to the integral equation based approximation (dashed red curve); b) free boundary 
positions computed for various mesh sizes; c) a solution profile H(x, r) for r = (blue line), 
t = T/2 (red curve), r = T (green curve); d) 3D plot and e) contour plot of the function H(x, r). 
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Table 2: Experimental order of convergence of the iterative algorithm for approximating the free 
boundary position. 



h 




eoc(L°°) 


err(i 2 ) 


eoc(L 2 ) 


0.03 


0.5 




0.808 




0.012 


0.215 


0.92 


0.227 


1.39 


0.006 


0.111 


0.96 


0.0836 


1.44 


0.004 


0.0747 


0.97 


0.0462 


1.46 


0.003 


0.0563 


0.98 


0.0303 


1.47 


0.0024 


0.0452 


0.98 


0.0218 


1.48 


0.002 


0.0378 


0.98 


0.0166 


1.48 
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T 

Figure 10: A comparison of the free boundary function q r (t) computed for the Risk Adjusted Pric- 
ing Methodology model. Dashed red curve represents a solution corresponding to R — 0, whereas 
the green curves represent a solution q r {t) for different values of the risk premium coefficients 
R = 5, 15, 40, 70, 100. 

4.2.2 Risk Adjusted Pricing Methodology model 

In the next example we computed the position of the free boundary g(r) in the case of the 
Risk Adjusted Pricing Methodology model - a nonlinear Black-Scholes type model derived 
by Jandacka and Sevcovic in [32]. In this model the volatility a is a nonlinear function of 
the asset price S and the second derivative dgV of the option price, and it is given by 
formula (O. In Fig. [10] we present results of numerical approximation of the free boundary 
position g R (r) = S R (T — r) in the case when the coefficient of transaction costs C = 0.01 
is fixed and the risk premium measure R varies from R = 5, 15, 40, 70, up to R = 100. We 
compare the position of the free boundary g R (r) to the case when there are no transaction 
costs and no risk from volatile portfolio, i.e. we compare it with the free boundary position 
q°(t) for the linear Black-Scholes equation (see Fig. [T0l . An increase in the risk premium 
coefficient R resulted in an increase of the free boundary position as it can be expected. 

In Tab.|3]and Fig.[l0]we summarize results of comparison of the free boundary position 
g R for various values of the risk premium coefficient to the reference position g = g° 
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Table 3 : Distance 1 1 q r ■ — g° \ \ p (p = 2 , oo) of the free boundary position g R from the reference free 
boundary position g° and experimental orders ctoo and a 2 of convergence. 



R 


IK -floe 




\\Q K -Q u h 


a 2 


1 


0.0601 




0.0241 




2 


0.0754 


0.33 


0.0303 


0.328 


5 


0.102 


0.33 


0.0408 


0.326 


10 


0.128 


0.33 


0.0511 


0.324 


15 


0.145 


0.32 


0.0582 


0.323 


20 


0.16 


0.32 


0.0639 


0.322 


30 


0.182 


0.32 


0.0727 


0.321 


40 


0.2 


0.32 


0.0798 


0.32 


50 


0.214 


0.32 


0.0856 


0.319 


60 


0.227 


0.32 


0.0907 


0.318 


70 


0.239 


0.32 


0.0953 


0.317 


80 


0.249 


0.32 


0.0994 


0.317 


90 


0.259 


0.32 


0.103 


0.316 


100 


0.268 


0.32 


0.107 


0.316 




20 40 60 80 100 20 40 60 80 100 

R R 



Figure 11: Dependence of the norms \\g R — g°\\ p (p = oo, 2) of the deviation of the free boundary 
g = g R (r) for the RAPM model on the risk premium coefficient B. 
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Figure 12: A comparison of the free boundary function q(t) computed for the Barles and 
Soner model. Dashed red curve represents a solution corresponding to R = 0, whereas the 
green curves represents a solution q(t) for different values of the risk aversion coefficient a = 
0.01,0.07,0.13,0.25,0.35. 




Figure 13: Dependence of the norms || g a — g°\\ p (p = oo, 2) of the deviation of the free boundary 
g = g a (r) for the Barles-Soner model on the risk aversion parameter a. 



computed from the Black-Scholes model with a constant volatility a = a, i.e. R = 0. The 
experimental order a p of the distance function \\g R — g°\\ p = 0(R ap ) has been computed 
for p = 2, oo, as follows: 

_ ln(||^- g || p )-ln(|| g ^- g °|| p ) 
ap In Ri - InRi-i 

According to the values presented in Tab. [3] it turns out that a reasonable conjecture on the 
order of convergence is that \\q r -q°\\ p = 0(i? 1/3 ) for both norms p = 2 and j> = oo. Since 
the transaction cost coefficient C and risk premium measure R enter the expression for the 
RAPM volatility © only in the product C 2 R we can conjecture that \\g ' — Q 0,0 \\ P = 
0{C 2 ^R^ 3 ) as either C -> 0+ or R -» 0+ 



36 



Daniel Sevcovic 



Table 4: Distance \\g a — g Q \\ p (p = 2, oo) of the free boundary position g a from the reference free 
boundary position g° and experimental orders oioo and of convergence. 



a 


II a Oil 

\\Q - Q loo 


Otoo 


II a U II 

\\g - g h 


OL1 


0.01 


0.156 




0.0615 




0.02 


0.25 


0.68 


0.0985 


0.68 


0.05 


0.472 


0.69 


0.184 


0.679 


0.07 


0.602 


0.72 


0.232 


0.69 


0.1 


0.793 


0.77 


0.298 


0.712 


0.11 


0.857 


0.82 


0.32 


0.74 


0.13 


0.99 


0.86 


0.364 


0.766 


0.15 


1.13 


0.92 


0.409 


0.807 


0.2 


1.52 


1. 


0.529 


0.897 


0.25 


1.97 


1.2 


0.669 


1.05 


0.3 


2.49 


1.3 


0.833 


1.21 


0.35 


3.07 


1.4 


1.03 


1.35 



4.2.3 Barles and Soner model 

Our next example is devoted to the nonlinear Black-Scholes model due to Barles and Soner 
(see [8]). In this model the volatility is given by equation (©. Numerical results are depicted 
in Fig. [321 Choosing a larger value of the risk aversion coefficient a > resulted in increase 
of the free boundary position £ a (r). The position of the early exercise boundary £ a (r) has 
considerably increased in comparison to the linear Black-Scholes equation with constant 
volatility a = a. In contrast to the case of constant volatility as well as the RAPM model, 
there is, at least a numerical evidence (see Fig{l2]and g a for the largest value a = 0.35) that 
the free boundary profile g a (r) need not be necessarily convex. Recall that that convexity 
of the free boundary profile has been proved analytically by Ekstrom et al. and Chen et al. 
in a recent papers [13, 19,20] in the case of a American put option and constant volatility 
a = a. 

Similarly as in the previous model we also investigated the dependence of the free 
boundary position g = g a (r) on the risk aversion parameter a > 0. In Tab. [4] and Fig. [13] 
we present results of comparison of the free boundary position g a for various values of the 
risk aversion coefficient a to the reference position g = g°. Inspecting values a p of the 
order of distance ||^ a — g°\\ p it can be conjectured that ||^ a — g°\\ p = 0(a 2 / 3 ) as a — > for 
both norms p = 2 and p = oo. 

5 Transformation methods for Asian call options 

Path dependent options are options whose pay-off diagram depends on the path history 
of the underlying asset. Among path dependent options Asian options plays an important 
role as they are quite common in currency and commodity markets like e.g. oil industry 
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(cf. [15,28]). Asian options may depend on the averaged path history in several ways. We 
shall restrict our attention to the so-called floating strike Asian call options. The floating 
strike price is assumed to be an arithmetic average of the underlying asset prices over the 
entire time interval [0, T] where T is the expiration time. 

Let us define the arithmetic average A = A t of the underlying asset S = St by 

1 I* 
A t = - S T dr. 

1 Jo 

For the case of the so-called Asian floating strike call option the pay-off diagram at expiry T 
reads as follows: V(S, A, T) = max(S' — A, 0). It means the price V of an option contract 
will depend not only on the underlying asset price S, time t, but also on the underlying 
asset path average A, i.e. V = V(S, A, t). 



5.1 Governing equations for Asian options 

As it is usual in the option pricing theory, we shall describe the asset price dynamics by a 
geometric Brownian with drift g, dividend yield q > and volatility a, i.e. dS = (g — 
q)Sdt + aSdW where W is the standard Wiener process. If we apply Ito's formula to the 
function V = V(S,A,t) we obtain 

„, (dV a 2 n 2 d 2 V\ dV 1C1 dV , . 

dV ={-m + Y S ds^) + ds dS+ dA dA - 

In the case of arithmetic averaging we have dA = t~ 1 {S — A)dt. Hence the differential dA 
is of the order of dt and this is why the above expression for dV indeed represents its lowest 
order approximation when taking into account stochastic character of the dynamics of the 
asset price S. Therefore, following standard arguments from the Black-Scholes theory one 
can derive the governing equation for pricing Asian option with arithmetic averaging in the 
form: 

dv <j 2 a2 d 2 v . .av s-Adv n 

Ht + Y S V + 5(r " q) OS + — 9A - rV = °' (75) 
where 0<t<T, S, A > (see e.g. [15]). For Asian call option the above equation is 
subject to the terminal pay-off condition 

V(S, A, T) = max(5 - A, 0), S, A > . (76) 

It is well known (see e.g. [15, 38]) that for Asian options with floating strike we can achieve 
dimension reduction by introducing the following similarity variable: 

x = j, W(x,r) = ~V(S,A,t) 

where r = T - t. It is straightforward to verify that V(S, A, t) = W(S/A, T - t)A is a 
solution of (1751 iff W = W(x, r) is a solution to the following parabolic PDE: 

dW a 2 2 d 2 W . . dW x-1 f TTT dW\ 



dr 2 



2 d 2 W , . dW x-1 f TTr dW\ 

*sr-<r-<)**r-T=?{ w -*&r) +rW - - <77) 
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where x > and < r < T. The initial condition for W immediately follows from the 
terminal pay-off diagram for the call option, 

W(x,0) =max(x- 1,0). (78) 
5.2 American style of Asian options 

Following Dai and Kwok [15], American style of Asian options is characterized by the 
exercise region 

£ = {(S, A, t) G [0, oo) x [0, oo) x [0, T), V(S, A, t) = V{S, A, T)}. 

In the case of a call option this region can be described by an early exercise boundary 
function Sf = Sf(A, t) such that 

£ = {(5, A, t) e [0, oo) x [0, oo) x [0, T), S > S f (A, t)}. 

For American style of an Asian call option we have to impose a homogeneous Dirichlet 
boundary condition V(0, A,t) = at S = 0. According to [15] the C 1 continuity condition 
at the point (St(A, t),A, t) of a contact of a solution V with its pay-off diagram implies 
the following boundary condition at the free boundary position Sf(A,t): 

^L(S f (A,t),A,t) = l, V(Sf(A,t),A,t)=S f (A,t)-A (79) 

for any A > and < t < T. It is important to emphasize that the free boundary function 
Sf can be also reduced to a function of one variable by introducing a new state function 
Xf(t) as follows: 

Sf(A,t) = Axfit). 

The function xf = Xf(t) is a free boundary function for the transformed state variable 
x = S/A. For American style of Asian call options the spatial domain for the reduced 
equation (|77T ) is given by 

0<x<q(t), tg(0,T), where q(t) = x f (T- r) . 

Taking into account boundary conditions (l79l for the option price V we end up with corre- 
sponding boundary conditions for the function W: 

dW 

W(0,r) = 0, W(x,t) =x-l, — (x,t) = 1 at x = q(t) (80) 

ox 

for any < r < T and the initial condition 

W(x,0) = max(x - 1,0) (81) 

for any x > 0. 
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5.3 Fixed domain transformation for American style of Asian call options 

Similarly as in Section 3, in order to apply the fixed domain transformation for the free 
boundary problem (1771 ), d80l ), (18TI ) we introduce a new variable £ and an auxiliary function 
II = n(£, r) (again representing a synthetic portfolio) defined as follows: 

e = ln^V m,T) = W(x,r)-x™. (82) 

Clearly, x G (0, q(t)) iff £ G (0, oo) for t G (0, T). The value £ = oo of the transformed 
variable corresponds to the value x = (S = 0) expressed in the original variable. On 
the other hand, the value £ = corresponds to the free boundary position x = q(t) (S = 
AS f (A,t)). 

A straightforward calculation similar to that of Section 4 enables us to to justify that the 
function n = II(£, r) is a solution to the following parabolic PDE 

an _ ,dn a 2 d 2 n ( 1 \ 

^ + ^ r >Bi-YW + { r + —r) n = - m 

where the term a(£, r) is given by 

, o(t) a 2 pe~t — 1 
^=^K + r- <7-- - ^= • (84) 

p(t) 2 r — r 

Notice the spatial dependence of the coefficient a = a(£, r) in comparison to the case of 
transformed equations for a linear or nonlinear Black-Scholes equations (see (l32l ). (l56l>). 
The initial condition for the solution II follows from (|8TT ) 

n(e ' 0) - \ £>lnp(0). 

Since W(£, r) = 1 for x = q(t) and W(0, r) = we conclude the Dirichlet boundary 
conditions for the function II 

n(0,r) = -l, II(oo,t)=0. 

It remains to determine an algebraic constraint between the free boundary function q(t) 
and the solution II. Similarly as in the case of a linear or nonlinear Black-Scholes equation 
we obtain, by differentiation the condition W(p(t),t) = p{r) — 1 with respect to r, the 
following identity: 

d . . dW . , . . d . . cW, , . . 

-j-py T ) = ~dx~^ T ' ,T ^ p ^ + ~q^^ t ^ t >- 

Since ^-(p(t),t) = 1 we have ^-(x,t) = at x = q{t). Assuming continuity of the 
function II and its derivative IT^ up to the boundary £ = we have 

2 d 2 w . . an . . dW . . 
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as x —* p(r). Passing to the limit x — > p(r) in equation (1771) we end up with the equation 



<r - q)p{r) - y^(0, r) + + r[p(r) - 1] = 0. 



It yields an algebraic nonlocal expression for the free boundary position q{t) 



= ' , , ' — • (85) 



+ T 



Next we determine the starting point of the free boundary function g(0). It means that we 
have to find the terminal value of the original state variable Xf(T). We shall assume a 
structural assumption 

r > q > (86) 

on the interest and dividend rates r, q. If g(0) > 1 then In g(0) > and this is why the initial 
function n(£, 0) is equal to —1 in some right neighborhood of £ = 0. Thus 9fn(£, 0) = 
for < £ <C 1. Again, assuming continuity of c^n(£, r) at (£, r) = (0, 0) we obtain 

an. , , sn, N on, . 

hm — -(0,r) = lim -^r(x,T) = hm — -(x,Q)=0. 



As a consequence we can conclude 



r + m 1 + rT 
g + i 1 + qT 



This initial condition for g(0) is exactly the same as the one derived recently by Dai and 
Kwok in [15]. They proved, for a general choice of r, q > that the initial condition for the 
function g is given by 

fl+rT 

0(0) = max — , 1 

^ V ' \l + qT 



In summary, we have transformed the free boundary problem for pricing American style of 
Asian call option with floating strike price into the following nonlocal parabolic PDE with 
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an algebraic constraint 

an _ N an a 2 d 2 u ( 1 \ „ 

/ \ 1+^-^ + ^-^ (0,t) ^ rp 

p{T) = i + g (r-r) ' 0<r<T ' 

subject to initial and boundary conditions 

n(0,r)=-l, n(oo,r)=0, 

n( ^ 0) -\ £>H(l + rT)/(l + q T)), (87) 
e (0) = (l + rr)/(l + gT), 

where 

5.4 An approximation scheme for pricing American style of Asian options 

Similarly as in the case of a nonlinear Black-Scholes equation (see Section 4) we restrict 
the spatial domain £ G (0, oo) to a finite interval of values £ G (0, L) where L > is 
sufficiently large. Let k > denote by the time step, k = T/m, and, by /i > the spatial 
step, h = L/n where m,n G iV again stand for the number of time and space discretization 
steps, resp. We denote by Hj an approximation of II(£j, Tj), p> w £?(r.,) where £j = i/i and 
Tj = jk. Then for the Euler backward in time finite difference approximation of equation 
(l60l we have 



IP -IP" 1 ,dIF /a 2 ^ e -*-l\dIP a 2 d 2 W ( 1 . 

+ - - + f-= ~ ir^rr- + \ r + - IP = 



k d£ \ 2 T - tj J d£ 2 d 2 ^ \ T 

(88) 

where iP is an approximation of the value b(jj) where the function 6(r) is defined as in 
d56l >, i.e. b(r) = + r — q. The solution IF = IP (a;) is subject to Dirichlet boundary 

conditions at £ = and £ = L. We set n°(f ) = II(£, 0) (see dSTJl. A gain we split the above 
problem into a convection part and a diffusive part by introducing an auxiliary intermediate 
step IP" 2 : 



(Convective part) 



(Diffusive part) 



Ip-2 - IP 



-1 



+ ^d x .rp-2=0, (89) 



IP'-IP'-I (a 2 p>e-t-l\dIP a 2 d 2 W . ■>,,,. 

Sr + ^ ^-^r^^+ r + IP = 0. (90) 



fc V 2 T-Tj J d£ 2 d 2 £ 

The convective part can be approximated by an explicit solution to the transport equation 
<9 r II + 6(r)c^n = for £ > and r G (r J -_i,r 7 ] subject to the boundary condition 
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II(0,r) = —1 and the initial condition II(£,Tj_i) = IF In contrast to classical 

plain vanilla call or put options the free boundary function q(t) need not be monotonically 
increasing (see e.g. [15] or [28]). Therefore depending on whether the term 6(r) = |© + 

r — q is positive or negative the boundary condition 11(0, r) = — 1 at £ = is either in- 
flowing (b > 0) or out-flowing (b < 0). It means that the boundary condition 11(0, r) = — 1 
can be prescribed only if 6(r) > 0. Let us denote by B(t) the primitive function to b(r), 
i.e. B(t) = In q(t) + (r — q)r. Solving the equation <9 r II + 6(r)3^Il = we obtain: 
ri(£,r) = IP'- 1 ^ - B(t) + Bfa-i)) if i - B(t) + Bfo-i) > and n(£,r) = -1 

otherwise. Hence the full time-space approximation of the half-step solution 2 can be 
obtained from the formula 

IT 7 " 2 " = / nJ_1 ^) tf ^ = Ci-ln^ +ln^_i- (r-q)k > 0, 
* \ — 1 otherwise. 

In order to compute the value IP -1 (77$) we make use of a linear approximation between 
discrete values II? -1 , i = 0, 1, n. 

Using central finite differences for approximation of the derivative d x IF we can ap- 
proximate the diffusive part of a solution d90~l ) as follows: 

II' I!' / 1 \ , 

~^k^~ + ( r + —J u > (92) 



2 T-Tj J 2h 2 /i 2 

Therefore the vector of discrete values IF = {II?, i = 1,2, ...,n} at the time level j G 
{1, 2, m} is a solution to a tridiagonal system of equations 

4^-1 , • '/I"; • "/111, II- (93) 



for i = 1,2, n, where 

oP. = a ^) = ___^ + __(_ + 



, , k 9 k ( o 1 p>e & — 1 



, , . 2 k ( a 1 p?e^ 1 - 1 , 



#' ee ^(^) = 1+ ^ r+ _L_^A ; -(Qj+^). 
The initial and boundary conditions at r = and x = 0, L, can be approximated as follows: 



n 



o 



-1 for& <]n((l + rT)/(l + qT)) 
for&>ln((l + rT)/(l + gT)) 
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for % = 0, 1, n, and Hq = — E, IFn = 0. The equation for the free boundary position 
g can be approximated by means of a finite difference approximation of d x TI at the origin 
£ = as follows: 

J = l + r(T- Tj ) + (T-^^ 

* l + q(T- Tj ) • ^ 

We formally rewrite discrete equations d9~T| ), d93l ) and d95l ) in the operator form: 

Q> = -F(IF), 

W-h = T(U j ,g>), (96) 

«4(^')ip = if'-s, 

where .F(IP) is the right-hand side of the algebraic equation d95l ), T(IP, £ J ) is the transport 
equation solver given by the right4iand side of (I9TI ) and A = A{q 3 ) is a tridiagonal matrix 
with coefficients given by d94l) . The system d96l > can be approximately solved by means of 
successive iterates procedure. We define, for j > 1, IP' = IP , g 3,0 = g 3 . Then the 
(p + l)-th approximation of IT- 7 and g> is obtained as a solution to the system: 

pj,P+i = jF(n^), 
n i-ip+! =T(n^)^> +1 ), (97) 
Aig 3 ^ 1 )!!^^ 1 = Tp-hp+ l . 

Now, if the sequence of approximate discretized solutions {(rP' p , g 3 ' p )} c ^ = i converges to 
some limiting value (IF' 00 , q>'°°) as p —> oo then this limit is a solution to a nonlinear sys- 
tem of equations (l96l > at the time level j and we can proceed by computing the approximate 
solution in the next time level j + 1. 



5.5 Computational examples of the free boundary approximation 

We end this section with several computational examples documenting the capability of the 
new method for valuing early exercise boundary for American style of Asian call options 
with arithmetically averaged floating strike. 

In Fig. [14] we show the behavior of the early exercise boundary function g(r) and the 
function Xf(t) = g(T — t). In these numerical experiments we chose r = 0.06, q = 
0.04, a = 0.2 and very long expiration time T = 50 years. These parameters correspond 
to the example presented in the preprint version of the paper [15] by Dai and Kwok. As 
far as other numerical parameters are concerned, we chose the mesh of n = 100 spatial 
grid points and we considered the number of time steps m = 2 x 10 6 in other to achieve 
very fine time stepping corresponding to 13 minutes between consecutive time steps when 
expressed in the original time scale of the problem. In order to make a graphical comparison 
of the state variable x = S/A to its reciprocal value A/S considered in [15] we also plot 
the function l/xf(t). 
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Figure 15: A 3D plot (left) and contour plot (right) of the function II (£, r). 
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Figure 17: A comparison of the free boundary position Xf(t) = g(T — t) (left) and l/xf(t) = 
1/ ' q(T — t) (right) obtained by our method (solid curve) and that of the projected successive over 
relaxation algorithm from [15] (dashed curve). 



In Fig. [I5]we can see the behavior of the transformed function II in both 3D as well as 
contour plot perspectives. Fig. [16] depicts the initial condition n(£,0) and five time steps 
of the function £ ^ IT(£, Tj) for tj = 0.1, 1, 5, 25, 50. A comparison of the free boundary 
position Xf(t) = g(T — t) as well as of its reciprocal value l/xf(t) = 1/ g{T — t) obtained 
by our method (solid curve) and that of the projected successive over relaxation algorithm 
from [15] (dashed curve) is shown in Fig. [17] It is clear that our method and that of [15] give 
almost the same result in the one third of the time interval [0, T] close to the expiration time 
T = 50. On the other hand, the long time behavior when r — > T is quite different. Notice 
that lim r _*x q(t) can be analytically computed and is equal to 1 (see [28] or [48]). Hence 
the long time behavior of g seems to be better approximated by our method. A comparison 
of early exercise profiles with respect to varying dividend rate q is shown in Fig. IT81 

Finally, we present numerical experiments for shorter expiration times T = 0.5833 
(seven months) and T = 1 (one year) with zero dividend rate q = and r = 0.05, a = 0.2. 
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Figure 18: A comparison of the free boundary position x / (t) = g(T — t) for various dividend yield 
rates q = 0.04,0.035,0.3,0.25. 




Figure 19: The free boundary position Xf(t) = g(T — t) for Asian call on assets paying no 
dividends (q = 0) with expiration time T = 0.7 (left) and T = 1 (right). 
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Conclusion 

In this survey paper we presented recent developments in fixed domain transformation 
methods applied to evaluation of the early exercise boundary for American style of op- 
tions. We discussed an iterative numerical scheme for approximating of the early exercise 
boundary for a class of Black-Scholes equations with a volatility which may depended on 
the asset price as well as the second derivative of the option price. The method consisted of 
transformation the free boundary problem for the early exercise boundary position into a so- 
lution of a nonlinear parabolic equation and a nonlinear algebraic constraint equation. The 
transformed problem has been solved by means of operator splitting iterative technique. We 
also presented results of numerical approximation of the free boundary for several nonlin- 
ear Black-Scholes equation including, in particular, Barles and Soner model and the Risk 
adjusted pricing methodology model. The method of fixed domain transformation has been 
also applied for evaluation of early exercise boundary for American style of Asian option 
with arithmetically averaged strike price. 
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